#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import opengate as gate
import opengate.geometry.materials as gmat
from pathlib import Path

if __name__ == "__main__":
    sim = gate.Simulation()

    # options
    # sim.visu = True
    sim.visu_type = "vrml"
    sim.random_seed = "auto"
    sim.number_of_threads = 2
    sim.progress_bar = True
    sim.output_dir = "./output"
    data_path = Path("./data")

    # units
    m = gate.g4_units.m
    mm = gate.g4_units.mm
    keV = gate.g4_units.keV
    Bq = gate.g4_units.Bq
    gcm3 = gate.g4_units.g_cm3

    # world
    world = sim.world
    world.size = [1 * m, 1 * m, 1 * m]
    world.material = "G4_AIR"

    # CT image
    ct = sim.add_volume("Image", "ct")
    ct.image = data_path / "ct1-4mm.nii.gz"
    ct.material = "G4_AIR"  # material used by default
    tol = 0.2 * gcm3
    mat_table = data_path / "Schneider2000MaterialsTable.txt"
    density_table = data_path / "Schneider2000DensitiesTable.txt"
    (
        ct.voxel_materials,
        materials,
    ) = gmat.HounsfieldUnit_to_material(sim, tol, mat_table, density_table)
    print(f"Reading CT image {ct.image}")
    print(f'Density tolerance = {gmat.g4_best_unit(tol, "Volumic Mass")}')
    print(f"Number of materials in the CT : {len(ct.voxel_materials)} materials")
    ct.dump_label_image = Path(sim.output_dir) / "labels.mhd"

    # Activity source from an image
    source = sim.add_source("VoxelsSource", "vox")
    source.attached_to = ct.name
    source.particle = "ion"
    source.ion.Z = 71
    source.ion.A = 177
    source.activity = 1e5 * Bq / sim.number_of_threads
    source.image = data_path / "spect-jm.mhd"
    source.direction.type = "iso"
    source.energy.mono = 0 * keV
    # compute the translation to align the source with CT
    # (considering they are in the same physical space)
    source.position.translation = gate.image.get_translation_between_images_center(
        ct.image, source.image
    )
    print(f"Reading source image {source.image}")

    # physics
    sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3"
    sim.physics_manager.enable_decay = True
    sim.physics_manager.set_production_cut("world", "all", 1 * m)
    sim.physics_manager.set_production_cut("ct", "all", 1 * mm)

    # add dose actor (get the same size as the source)
    source_info = gate.image.read_image_info(source.image)
    dose = sim.add_actor("DoseActor", "dose")
    dose.output_filename = "ctdose.mhd"
    dose.attached_to = ct.name
    dose.size = source_info.size
    dose.spacing = source_info.spacing
    # translate the dose the same way as the source
    dose.translation = source.position.translation
    # set the origin of the dose like the source
    dose.output_coordinate_system = "attached_to_image"
    dose.hit_type = "random"
    dose.dose_uncertainty.active = False
    dose.dose.active = True

    # add stat actor
    stats = sim.add_actor("SimulationStatisticsActor", "Stats")
    stats.track_types_flag = True
    stats.output_filename = "stats.txt"

    # go
    sim.run()

    # end
    print(f"Output statistics are in {stats.get_output_path()}")
    print(f"Output edep map is in {dose.edep.get_output_path()}")
    print(f"vv {ct.image} --fusion {dose.edep.get_output_path()}")
    print(stats)
